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Second-order and fourth-order two-dimensional advective equations 
developed by W. P. Crowley were examined to determine their applica- 
bility to atmospheric models. 

One second-order and two fourth-order forms were evaluated from 
their performances on a simple pattern advec.ted.-hy_ a linear, divergent 
velocity field. The same equations were substituted for the advection 
term of a simple barotropic forecast model to determine their perform- 
ances on more general non-linear conditions. 

All forms of the equations remained stable in time and demonstrated 
the phase and amplitude characteristics predicted by Crowley. The 
fourth-order "advection" form gave best results. 

When substituted in the barotropic model, the fourth-order forms 
lead to improved trough and low-center movements, but RMSE was slightly 
larger than that resulting from second-order forms. The better RMSE 
of the lower order forms apparently resulted from their diffusive 
charateristics. 
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CHAPTER I 



INTRODUCTION 



The hydrodyamical equations governing fluid motion usually contain 
a non-linear advective or transport term which often plays an important 
role with respect to the temporal changes of the pertinent parameters. 
It is, therefore, desirable to approximate the advective term accu- 
rately when a numerical solution of the equations is sought. 

In order to investigate the accuracy of various finite-difference 
representations of the transport term, it will be assumed that some 
parameter, <j>, is conserved in a two-dimensional flow field. This 
assumption may be expressed in Eulerian form as 



H = - V • V«j>, (1.1) 

-V 

where V is the horizontal velocity and v is the two-dimensional del 
operator* 

When written in scalar form, eq. 1.1 becomes 



a£ 

8t 



,.3i +u M 

“ 3x + V ty ] ’ 



( 1 . 2 ) 



where u and v are the x and y components, respectively, of the vector 
velocity. If the flow is non-divergent, eq. 1.2 may be expressed as 



|| - - 0(¥,<l>), (1.3) 

where n is the stream function and J is the well-known Jacobian 
operator. 

An expression equivalent to eq. 1.2 is 

M. s faCcpu) . a(<t>v) ,'au av n 

at ' ax ay " * ax av * U * 4J 
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Although the analytical solution to eq. 1.4 is identical to that of 
eq. 1.2, the latter form leads to a different finite-difference expres- 
sion possessing certain advantages which are discussed in Chapter II. 

For convenience, the advective equation as expressed by eq. 1.2 
will be called the "advection" form; equation 1.4 will be referred to 
as the "conservation" form. 

Although the analytical solution of a differential equation is 
exact, the finite-difference representation is only an approximation. 
When modern digital computers are utilized, the accuracy of the 
approximation is principally dependent on (1) truncation error, and 
(2) computational stability. Truncation error may be reduced (but 
never eliminated) by retaining high-order terms and by reducing the 
size of the grid-mesh. Elimination of instability error is more 
difficult because computational stability is dependent on the physical 
properties of the fluid as well as on the specific form of the finite- 
difference equation. 

Many finite-difference techniques have resulted from attempts to 
control instability while reducing truncation error. The purpose of 
this study was to evaluate by numerical experimentation the appli- 
cability to the atmosphere of high-order, two-dimensional advective 
techniques proposed by Crowley [1]. 

Examples of practical problems which require accurate advective 
techniques are (1) radioactive debris tracking, (2) air and water 
pollution studies, (3) estimation of movement and concentration of 
potentially lethal gases, and (4) numerical weather prediction 
Numerical models which provide solutions to these problems often 
require a high degree of accuracy throughout long periods of time. 
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For example, radioactive debris projected into the stratosphere may 
travel around the earth several times before the combined effects of 
turbulent diffusion, velocity divergence and deformation, and fall-out 
reduce the concentration of the radioactive substance to an insignifi- 
cant level . 

A good advection scheme, therefore, must exhibit not only small 
errors in each time step, but also small cumulative errors after many 
steps. To determine the extent to which Crowley's equations satisfy 
this requirement, the initial experiment was designed to evaluate the 
performance of the techniques, given a known configuration of a 
conservative quantity and a known velocity field.- From observations 
of the position, shape, and magnitude of ' the quantity after many hours 
of advection, statistics on effective phase speed, conservation of 
intensity and shape, and computational stability were obtained. 

Since diagnostic and prognostic models of the atmosphere must 
approximate reaf conditions in a satisfactory manner, the advective 
terms employed must demonstrate stability and accuracy under non- 
linear conditions. The characteristics of Crowley's equations, given 
non-linear conditions, were evaluated by comparing the results obtained 
from a simple barotropic forecast model in which appropriate forms of 
the various advective equations were utilized- 
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CHAPTER II 



THE ADVECTIVE EQUATIONS 

The alternate forms of the advective equation will be developed in 
this chapter. Both result in a forward time step based upon centered 
space differences. Since the resulting finite-difference equations 
will be used eventually in the barotropic model, absolute vorticity, n, 
will be substituted for <j> in eqs. 1.1 - 1.4. 

The following is an explanation of the notation used to indicate 
time and space differences: 

1. Superscript n: denotes the current hour. 

2. Subscript i: denotes the x-coordinate of the grid point, 

x = iAX. 

3. Subscript j: denotes the y-coordinate of the grid point, 

y = JAy. 

The development of equations 1.2 and 1.4 follows closely that of 
Crowley with minor differences in notation. 



THE ADVECTION FORM 



A Taylor series expansion in time of equation 1.2 in one dimension 
results in 



n+1 

n 

i 



n 

0 + 

i 



At 



9n 

at 



n 



(At) 2 a 2 n 

2 at 2 



+ 0(At 3 ). 



( 2 . 1 ) 



From equation 1.2, 



In = n M 
at ax • 



( 2 . 2 ) 



Differentiating the above equation with respect to time gives 



a 2 n 


au an a 


[3nf 


at 2- = ' 


at ax at 


[ax] 
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If the order of differentiation is changed in the last term, and 



substitution is made from eq 2 2, the result is 



5' n .. 
3t* 



3u 3n . .. 3 
at ax 9x 



which may be rearranged to give 



3 2 r, 

3t 2 



,.2 3 2 n 

u dW “ 3x 



3n 



3u 

3t 



■ u Si 



u 3J1| 

3xJ 



(2.3) 



(2.4) 



Substitution of eq. 2.2 and 2,3 into 2,1 gives 

‘ >+1 n” ,,‘t H“. ( “ At)2 ^ " (it)2 8nl °f 311 U 3 “' ” 

n * 1 - uAt gj-, * — 2 3 xT r~ 3X| 3t U 3x 

i i ‘i i i 

If fx and dx^ are eva ' uated to second order and the final term of 
eq. 2.4 is discarded, substitution of the finite-difference forms 
results in 



n+1 



n 



n 1 


UAt] 




IjuAt 


2 ' 


n - 2 
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AX 1 


n -n 

i+l i-1' 


+ 2l AX 

V j 





n n \ 

n - 2n + n 
■1 i i-1' 



(2.5) 



which is second-order accurate in space and will be referred to as the 

"second-order" form Note, however, that if u is not constant, the 

equation is only first-order accurate in time. 

It is possible, although tedious, to develop, by similar methods, 

a form of eq. 2.5 which is fourth-order accurate in space. However, 

the "fourth-order" form so obtained remains only first-order in time. 

The finite-difference form is 
n+1 n , rnulf f n n 1 f n n 



n 



12! 


'Ax J 


_i| 


UAt] 


24 ! 


Iax J 
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f uAt] 


12 


AX 

V. J 


_1 


j UAt . 
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n r n 


n 1 


' n 
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On - 16 in 


+ n + 


ri 


+ n 


i 1 i+l i-P 


1 i+2 


i-2'_ 


f n n 


( n 


n 7 




|n - n 


- n - n 




i+l i-1' 


1 i+2 


±-2 } _ 





6n - 4 j n 



n 



i n ] + fn 11 + n n t 



2HAX] Q T’ i+1 i.J ( i+2 i- 2 JJ 
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( 2 . 6 ) 



THE CONSERVATION FORM 



The one-dimensional form of eq. 1,4 is 

dp 9u 3(nu) 

3t " n 8x ' 3x ( 2 / 7 ) 



The transport term, s(mi)/3x, can be considered to be a divergence of 

a flux. By Green's theorem, any change in the advected quantity with 

time is proportional to the differential of flux across the boundaries. 

If the boundaries are established at points half the distance between 

grid points, then the flux at one of the boundaries is expressed by 

x f i+% 

W uit 

If the variation of n is assumed to be linear within the zone, then 
integration of eq. 2.8 produces 




(2.9) 



If the variation of n is described by a cubic equation centered on the 
point, i + k, the resulting equation is 



P At - 1 fuAt] Fqf 
'i-t-% ax T6 Iax J v 



n i+l + n i^ " ^ n i+2 + n i-l^ 



I AX ] j^ n i+l" n i^ ” ^ n i+2" n i-l^ 



4 1 * 11 



l 

TZ 



uAt) 
AX I 



If UAt) 

+ W [ET, 



(n i+ i+ n ± ) - (n i+2 + n^) 

3 ( n i+l" " ^ n i+2" n i-l^ 



( 2 , 10 ) 



where the equation has been multiplied by and • is evaluated at 

i + H. 
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When Is evaluated to second-order and the flux divergence is 

dX 

substituted for the transport term, the second-order difference form 
of eq„ 2.7 is 



ni-i 



nr 

n L 1 

i L 



Hf (u i+r u i-i ) j' 

X T ^ 



n 



- F ) 
i-% 



( 2 . 11 ) 



where the form of the flux term is given by eq. 2,9. 

Evaluation of to fourth-order and substitution of flux diver- 

oX 

gence as computed by eq. 2 10, results in the fourth-order finite- 
difference conservation form 



n+1 

n 

i 




TZAx {8(u i+i' 



J i-l) " ^ U i+2” 




AX [ 




STABILITY CHARACTERISTICS 



( 2 . 12 ) 



It is convenient to express the advective equations in terms of 
linear operators. Both forms of the advective equations described in 
this chapter have tne same general equation, 



n-t-1 



1 n • i n ) 


n 


In j - An = 

{ J 1 l Ji 


(I-A)n 



(2.13) 



where I represents the identity matrix and A represents the particular 
linear difference operator. Stability is, therefore, dependent on the 

eigenvalues of (I - A). Crowley [1] has shown that eqs. 2.5, 2,6, and 

fuat) 2 



2.11 are stable for 
1.5. 



UAt 
< 

AX - 



1, and that eq. 212 is stable for 



AxJ - 



From Crowley's analysis of the amplitude and phase errors of the 
advection and conservation equations, the following generalizations 
are made: 
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1 Maximum amplitude damping and minimum phase retardation both 

occur at ^7 =0*7, 
ax 

2. The rate of decrease in damping and of increase in phase 

retardation is greater for values of larger than 0.7 than 
it is for values smaller than 0.7 
3 Amplitude damping and phase retardation are serious only for 
wave lengths shorter than four grid intervals. 

4, The conservation form has slightly greater errors than the 
advection form in both amplitude and phase. 

5. The fourth-order forms are clearly superior to the lower order 
forms in both amplitude and phase, 
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CHAPTER III 



EXPERIMENTAL PROCEDURES 



Both forms of the advective equation were developed for one- 
dimensional systems. Two-dimensional systems may be constructed by 
employing Marchuk's method of fractional time steps [1,3]. The linear 
operator form of such a system is 



for the conservation form. An and Bn represent the fluxes in the x 
and y directions respectively, and k represents the velocity diver- 
gence term in the conservation equation. 

Using Marchuk's method, the final forms of the advective equations 
may be expressed as follows: 

Advection form: 



time n + 1, which results from the first one-dimensional step. Since 
each fractional step is, by itself, stable, the combination of the steps 
retains stability. 




(3.1) 



for the advection form, and 



n+1 



n 



n 



* ( 1 -k) ( I-A) ( I-B)n 



(3.2) 



n* , - (I-B)n 

1 > J L 



n 




(3.3) 



Conservation form: 





(3.4) 



* X 

The symbol, n , represents an intermediate value of the quantity at 
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EQUATION CHARACTERISTICS TEST 



The purpose of the initial experiment was to determine quantita- 
tively, the performance characteristics of the several equations under 
known conditions which have some similarity to the atmosphere. 

The grid system employed consisted of 3,969 points arranged in a 
square matrix of dimensions 63 x 63. The interval, d, between grid- 
points was approximately two hundred nautical miles. In this grid, 
the point with coordinates (31,31) lies exactly at the center, so that 
distances from the center are proportional to [ ( 31 -i ) 2 + (31-j) 2 ]^. 

For simplicity in this experiment, the surface was assumed to be 
flat; therefore no map factor was introduced. 

A velocity field was generated which was proportional to the 
distance from the center point and which contained no shear (solid 
rotation). A pure divergent component was then added to the rotation. 
The resulting equations are 



where x = ( i - 3 1 ) and y - (j-31) are dimensionless coordinates. The 
first terms on the right are the rotational components and the second 
terms are the divergent components. M and N are arbitrary constants 
with units of meters per second, which determine the magnitude of the 
velocity components. 

The initial field of the quantity to be advected consisted of a 
sinusoidal perturbation centered on grid-point (31,19) and zeroes 
elsewhere The equation of the perturbation is 



u = - My + Nx , 



and 



v = Mx + Ny 



(3.5) 




(3.6) 
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where R = distance In grid intervals from point (31,19). Since d is 
approximately two hundred nautical miles, the perturbation has a 
maximum wave length of approximately two thousand nautical miles, 
which is of the same order of magnitude as the large scale vorticity 
disturbances . 

The computer program was designed so that the angular velocity and 
the divergence could be varied in sign and magnitude as well as the 
number and length of time steps. Values representative of the atmos- 
phere and forecast models were selected. The equation which relates 
the variables is 

27iR c d = MR c S(At), (3.7) 

where 

R c = distance from grid center in grid intervals = 12, 
d - length of grid interval = 3.81 x 10 5 meters, 

S = number of time steps. 

At = length of time step, 

M - velocity factor (meters per second). 

Preliminary calculations determined that with one-hour time steps, 
approximately two hundred steps would be required to advect the 
quantity through 2ir radians with a realistic wind. The necessary 
velocity factor is 



2nd - (6.28 )(3. 81x10 s ) • 
“ SAt " (2xl0 2 ) (3, 6x10 3 ) " 



3.32 m./sec. 



(3.8) 



This factor produces a linear velocity of approximately forty meters 
per second at the point (31,19). 

Initial experiments indicated that a suitable magnitude of N, the 
divergent factor, was 0,5 meters per second. 
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The scheme of the experiment was to advect the pattern through ir 
radians with positive divergence, then to change the sign of the 
divergence for the remaining it radians. Assuming perfect advective 
properties, this method would return the pattern to its initial posi- 
tion and configuration after advection through 2ir radians. Comparison 
of the initial and final patterns would show quantitatively errors in 
phase velocity and amplitude, as well as any computational instability 
which developed 

For purposes of comparison, eq. 1,3 was also subjected to the 
above conditions. The finite-difference form of eq. 1.3 which was 
used was 

(3.9) 

which is second order in both space and time. 



n+1 



n-l 



n i,j = n i,j 



mAt 

” 



u . 



n 
,-n. 



+ v. 



n 

n. 



n 



,-n, 



i.jv i-i-l.j i»j( i»j+l 



THE BAROTROPIC EXPERIMENT 

The experiment just described subjects the advective equations to 
velocities which are linear in form and in comparison to the atmos- 
phere, rather idealized. In order to test the advective equations 
under more realistic conditions, a simple form of the barotropic 
forecast model was chosen as follows: 



n2 H jjfi 9! _ c 

v at " at ' " V 

where 



(3.10) 



f = stream function 
y = constant = 4 

f = average value of coriolis parameter = 1.03126 x 10 -4 sec. -1 
Z = average value of 500-MB height = 5574 m. 

F g = advective forcing function = V * vn . 
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The form of the linear balance equation employed to compute the 

h f ' ' • ft 

stream function was 

V 2 Y + ^(w-Vf) = £ V 2 Z„ (3.11) 

However, the second term on the left in equation 3.11 was approximated 

by 

i(l VZ Vf), 
f f 

reducing eq 3 11 to a Poisson-type equation. 

The scheme of the experiment consisted of equating the various 
forms of the advective equations to F . This could not be done 

a 

directly since F. - V * V n, but by subtraction of Iq n from both sides 

a 

of eq. 2.13, the desired quantity was obtained as follows: 

n+1 n n 

q - In = - An = An. (3.12) 



The quantity Aq n is substituted for F in eq. 3.10. 

a 

Two different methods were used to integrate the stream function in 
time. When the Jacobian was used as the forcing function in the nrog- 
nostic equation, time integration of the stream function was accomplished 
by the leap-frog method. When the advection and conservation forms were 
used, time integration was done by forward steps, i.e., 



n+1 n 

y = y 



3t 



At 



(3.13) 



which is equivalent to eqs. 3.1, 32, 3.3, or 3.4. 

■i . ■ ■ 

PROGRAMMING PROCEDURES 

All numerical computation required by this study was accomplished 

' L T i _ * 

on the Control Data 6500 computer installed at Fleet Numerical Weather 
Central. FORTRAN, as modified for the CDC 6500, was the programming 
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language used. In addition to FORTRAN, several FNWC library sub- 
routines were utilized. 

No particular effort was expended to optimize the programs. 
Nevertheless, relative time requirements of the various methods should 
be significant. 

During the construction of the computer program to solve the 
finite-difference equations, several minor modifications of the 
advective equations were necessary. 

First of all, the boundaries of the grid necessitated the establish- 
ment of boundary conditions. The method chosen for the quantity, n, 
was simply to hold the boundary values constant. In the case of the 
fourth-order forms, the first interior value was computed using a 
second-order form. 

Since the flux value was evaluated at' half grid-points, and the 
boundary values of n were not allowed to' vary, the first computation 
of flux occurred between the boundary point and the first interior 
point. This was accomplished with a second-order form in the fourth- 
order equations. 

Similarly, the u and v components of the velocity field for the 
baro tropic experiment were computed by second-order methods near the 
boundary, but in a different manner. 

Since the u and v components are derived from the gradient of f, 
it was possible to use fourth-order methods along the boundaries which 
were perpendicular to the separate velocity components. For the 
boundaries which were parallel to the components, a second-order form 
was used for the first interior value. This particular method was 
necessary because the conservation form includes interpolated velocity 
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components and g divergence t§rm which requires the first derivatives 
qf the velocity components. 

The remaining modifications were due to the fact that the flux 
divergence term was computed for the points i + % and j +_ Since 
the y and V components were computed for the interger points only, it 
Wgs necessary to Interpolate tor the values at the intermediate points. 

This was accomplished by using a linear interpolation op the 
perpendicular boyndaries and a secondrorder interpolation elsewhere. 

fig, 1 shows the grid system employed ip tpe conservation method. 
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Figure 1. Schematic of grid used in the conservation 
method. Integer grid-points are represented 
by + ; half- integer points by o. 
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CHAPTER IV 



EXPE RIMENTAL RESULT S 

Ail contoured figures which appear in this chapter were produced 

by a California Computer Products, Inc , incremental plotter. The 

instructions to the plotter were computed by the CDC 6500 computer by 

a quadratic interpolation of grid-values Contour labels on figures 

applicable to the linear advection experiment are proportional to the 

actual values The contour origin, 0, corresponds to an actual value 

of -1.0 and the contour interval corresponds to 1.5 units 

To facilitate the discussion of the various advective equations, 

the following symbols will be used: 

C4 denotes four th-order conservation 
A4 denotes fourth-order advection 
A2 denotes second-order advection 
J1 denotes common Jacobian 

EQUATION CHARACTERISTICS TEST 

The experiment consisted of advection of a pattern through 2 tt 
radians, about the center grid- point This was accomplished in two 
hundred’ time steps of one hour each The initial pattern and the 
resulting patterns after advection through u and 2 rr radians were 
plotted Figs 2-4 show the results when the fourth-order conserva- 
tion equation was used. Fig. 5 and fig 6 show the final results when 
the advection equations were used. 

All methods show good results, considering that the final pattern 
represents' the cumulative errors in advection of the pattern once 
around the earth at a velocity representative of a moderate jet stream. 
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Figure 2. Contour plot of initial perturbation. Contour interval cor 
responds to 1.5 units; contour label 06 corresponds to 0.5 units. 
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Figure 3. Plot of perturbation after 100 time steps by C4 method. 
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Figure 4. Plot of perturbation after 200 time steps by C4 method. 
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Figure 5. Plot of perturbation after 200 time steps by A4 method. 
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Figure 6. Plot of perturbation after 200 time steps by A2 method. 
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There are, however, significant quantitative differences Table I 
shows the percentages of maximum error in phase and amplitude of each 
of the methods. 

For purposes of comparison, the experiment was repeated using the 
J1 equation Prior to the completion of one hundred steps, instability 
near the grid boundaries completely destroyed the pattern. This 
instability was eliminated by using one-half hour time steps. So 
that some reasonable comparison could be made between the methods, the 
number of time steps was doubled. The results are shown in fig. 7. 

A particularly interesting result of the latter test was that 
after advection over approximately equal distances using twice the 
number of time steps, the amplitude error of the J1 result was less 
than that of the A2 method On the other hand, the phase error was 
larger than that resulting from the A2 method 

All of the results were qualitatively consistent with Crowley's 
theoretical computations which were summarized in Chapter III. The 
eccentricity of the resultant patterns was apparently due to the wave- 
length sensitivity mentioned by Crowley. Support for this statement 
was provided by a comparison of the initial and final positions of 
the outermost contour in the forward semicircle; they were nearly 
coincident in both fourth-order results. 

Although it is not apparent from the figures, small -amplitude 
noise was created by the advection and conservation methods. Nowhere 
did the amplitude of this noise exceed five percent of the magnitude 
of the perturbation. Except for the area close "behind" the pattern, 
the magnitude of the noise was less than one-half percent that of the 
perturbation 
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TABLE I 



PERCENTAGE OF MAXIMUM ERROR IN PHASE 
AND AMPLITUDE OF 
ADVECTIVE EQUATIONS 



ADVECTIVE 

EQUATION A4 C4 A2 



PHASE 

ERROR 0.67% 1.99% 3.45% 



AMPLITUDE 

ERROR -25% -34% -51% 



J1 



4.38% 



-48% 
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Ftgure 7. 



Plot of perturbation after 400 one-half hour time steps 
by J1 method. 

37 



A variation of the basic experiment was conducted to determine 
whether the length of the time step would have an effect. The time 
step was reduced to one-half hour and the number of time steps was 
doubled. The results obtained using the C4 method were nearly identical 
to those of the basic experiment. In this case, at least, any effect of 
reduction of the time step was apparently counteracted by the increased 
number of steps. 

In his conclusions, Crowley [1] surmised that given a divergent 
velocity field, the performance of the conservation equation would be 
superior. He attributed the anticipated superiority to the explicit 
divergence term in the conservation form in contrast to the implicit 
divergence of the advection form. This experiment has not substan- 
tiated that supposition. 

Computation time is usually a consideration in the design of 
atmospheric models. Obviously, the high-order schemes treated in this 
study require more computer time than the mathematically more simple 
schemes. Due to the "purity" of the computation in this experiment, 
it is possible to estimate, with high confidence, the relative require- 
ments of the various forms by simply comparing the total computation 
time required. 

If the time required by the J1 method is considered to be unity, 

then the relative "cost" of the other methods is listed below. 

C4 = 5.14 
A4 = 3.74 
A2 = 1.98 

It is obvious that the cost of the high-order equations in a pure 
advective model is considerable, but when these equations are used in 
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an iterative model such as the barotropic forecast model, the percentage 
increase in total time is small 

THE BAR0TR0P1C EXPERIMENT 

The 500-MB analyses used In this experiment were produced by the 
Fleet Numerical Weather Central operational upper-air analysis 
program Eight consecutive analyses beginning with 0000Z 15 January 
1968 and ending with 1200Z 18 January were processed. 

The choice of period was random except for the requirement that 
it be a winter period. As it turned out, it may have been an 
unfortunate choice because the upper air codes were changed on 
1 January 1968. A comparison of the average number of 500-MB reports 
successfully decoded during the period of study with the average number 
successfully decoded during September 1968 results in the ratio, 

330/374, It is possible, therefore, that one or more of the analyses 
contained significant errors. 

Statistics on RMSE, pillow, variance of absolute vorticity, and 
mean absolute vorticity were computed for each map time. Contour maps 
of forecast values, and associated error maps, were produced for the 
following times only: 

1200Z 15 January 

0000Z 1 7 January 

1200Z 18 January 

The choice of these particular map times resulted from efforts to 

represent the extremes and mid-point of the period, and to avoid a 
multiple of twenty-four hours. 

In addition to the statistics based upon the results of the various 
advective equations, RMSE and pillow were computed from the operational 
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prognoses of FNWC. Contour plots of the FNWC operational forecasts 
corresponding to the three time periods listed above were also 
produced. 

Forecasts for twenty-four and forty-eight hours were produced. 

A one-hour time step was employed and no smoothing was applied to 
the forecasts produced by the advection and conservation equations. 

Early in the experiment, it was found to be necessary to smooth 
the stream fields when employing the 01 equation, due to short-wave 
instability near the boundary. The method adopted was to apply a 
heavy smoother to the first four interior rows of grid-points on the 
initial stream field and on the forecast stream field at the twenty- 
four hour point. The boundary points and other interior points were 
not smoothed. The equation of the smoother is 

i,j 4 1 ; 

This smoother was probably excessive, but no significant activity 

occurs south of 13°, the most northerly latitude of the smoothed points. 
The smoothing resulted in a significant reduction of c^ 2 , however. 

Table II shows the RMSE which resulted from the forecasts pro- 
duced by the C4, A4, A2, and J1 equations, and the FNWC operational 
forecasts. Means and standard deviations for the eight map times are 
also shown. 

Although RMSE is certainly not a' conclusive measure of forecast 
skill, it is usually an indication of the relative accuracy of various 
methods. Table II suggested that in this experiment, skill was 
inversely proportional to the complexity of the advection equation 
employed in the model. This was indeed disturbing. Furthermore, it 
is surprising to note that all methods except C4, without smoothing. 
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RMSE IN METERS OF VARIOUS FORECAST METHODS 
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surpassed the skill of the FNWC ODerational model, which is considerably 
more complex than the model utilized in this experiment. 

To resolve this apparent paradox, the various prognostic charts 
and their associated error charts were compared with the verifying 
analyses. From this analysis the following conclusions were drawn: 

1. Closed lows existing on the analysis continued to deepen 
when C4 and A4 equations were used. When all methods 
continued the deepening process the C4 and A4 equations 
showed the highest rate of deepening. Frequently, these 
lows showed filling on the verification maps. 

2. Movement of low centers was more rapid in the C4 and A4 
results, and usually more accurate when the low moved rapidly. 

3. Trough movements were significantly larger in the fourth-order 
results than they were in the results of the low-order equa- 
tions. Frequently, the trough movement forecast of the A4 and 
C4 equations exceeded the actual movement in the vicinity of 
the jet stream. 

4. None of the methods were successful in the development of a 
low-latitude cut-off low; a typical barotropic model failure. 

5. The southern extremity of all short-wave troughs, except those 
of FNWC, showed considerable lag. Considering this fact, the 
overall movement of troughs appeared to be the best in FNWC 
progs. Note, however, that FNWC has incorporated a modifica- 
tion tested by Lewit [4] to accelerate low-latitude troughs. 

From these conclusions, the inverse relationship of RMSE is at 
least partially explained. They do not, however, explain why RMSE of 
FNWC progs was large. Examination of the error charts showed that in 
general, FNWC error centers were smaller than those in the experimental 
methods. In the Pacific Ocean west of the Hawaiian Islands and over 
eastern Asia, however, FNWC errors were consistently larger than those 
in the experimental methods. A large positive error center in the 
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central Pacific and a moderate negative error over eastern Asia 
persisted on all FNWC charts „ The Pacific error center appeared on 
charts of the experimental methods, but was of lower intensity,, The 
negative error over Asia was absent in the experimental methods. No 
explanation was found for this inconsistency, 

Although generalizations were helpful in the explanation of the 
RMSE, they may not represent the characteristics of the model in a 
given forecast situation. A detailed analyses of the 15/1200Z map 
time was conducted to determine the qualitative and quantitative 
differences in performance of the models 

The 15/1200Z map time was chosen because (1) the report count was 
the largest of all map times in the series, (2) the report count for 
the verification time, 16/1200Z, was well above average, and (3) the 
RMSE of the various methods closely approximated the mean values for 
all eight periods. 

Fig, 8 is a contour plot of the 15/1200Z analysis. It is typical 
of a rather low-index pattern with a nearly symmetrical blocking situa- 
tion in the central Pacific and deep troughs over the eastern United 
States, eastern Atlantic, and eastern Europe 

Fig, 9 is a plot of the 16/1200Z analysis, which provided the 
verification for the prognoses. It is apparent that little change 
has occurred in the Pacific hemisphere, but that the movement of 
short waves in the Atlantic hemisphere has caused changes sufficient 
to exercise the forecast models. 

Figs. 10-14 are plots of the twenty-four hour progs produced by 
each of the methods. The most striking feature of these figures is 
the smoothness of the advection and conservation methods compared to 
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Figure 8. FNWC height analysis at 500-MB for 1200Z 15 January 1968. 
Contour interval is 60 meters. 
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Figure 9. FNWC height analysis at 500-MB for 1200Z 16 January 1968. 
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Figure 10. 



C4 24-hour 500-MB prog verifying 1200Z 16 January 1968. 
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Figure 12. A2 24-hour 500-MB prog verifying 1200Z 16 January 1968. 
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Figure 13. J1 24-hour 500-MB prog verifying 1200Z 16 January 1968. 
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the J1 method. This is a consistent feature. The smoothness of the 
FNWC prog is a result of repeated smoothing during the integration 
process. 

The trough over Labrador and associated ridges to the east and west 
appear to have been handled best by the C4 method, The FNWC prog is 
nearly as good in position and shape but has a much larger error in 
central pressure. This is a deepening situation. 

Large differences occurred in the handling of the low which, in the 
analysis, is located off the northeast Labrador coast. This low filled 
by verification time The J1 method and FNWC excelled in this case. 

Note that the Crowley methods moved this feature northeastward and the 
A4 and A2 methods deepened it Contrary to the generalization made 
earlier, the C4 method showed slight filling. 

The trough in the eastern Atlantic was handled rather well by all 
methods, but the J1 method appears to have the best shape. Note that 
the A4, C4 and FNWC methods all move the trough too rapidly in the 
higher latitudes. 

The weak low located over southern Norway is of particular interest, 
since it is the only rapidly moving low on the chart It can be fol- 
lowed to a position approximately 300 miles north of the Black Sea by 
verification time. The movement forecasted by the C4 method is 
clearly superior, but the central pressure is in error by more than 
120 meters, which is more than twice the error shown in the FNWC prog 

The deep trough near the Caspian Sea was handled poorly by all 
methods. The fourth-order methods are particularly bad because of 
excessive deepening, especially the C4 method. Central Asia was 
handled similarly by all methods with only small differences in detail. 
The FNWC prog, however, appears to be overly smooth 



51 



Although there appeared to be but small differences in the handling 
of the deep low north of Japan, the strong gradient produced quite 
different error patterns. The low moved northeast and deepened. 

The J1 , A4 and A2 methods moved it north. The FNWC prog resulted in 
movement slightly west of north. The C4 prog was clearly superior in 
the movement and deepening of this center, but also had an excessive 
northerly component. 

Since the western Pacific area was particularly difficult for the 
FNWC prog, it was interesting to compare the error charts of the FNWC 
prog with the C4 prog in this area. Fig. 15 and fig. 16 are computa- 
tions of prog minus verification, contoured at intervals of sixty 
meters. The differences are quite significant. 

All methods treated the blocking situation similarly. All failed 
to reduce the intensity of the low to the value shown on the verifica- 
tion. Note, however, that the low is positioned in a sparse data area. 
It is quite possible that the center was more intense than indicated 
by the analysis. 

No significant differences in the handling of the Gulf of Alaska 
cyclone were apparent. 

The foregoing analysis of movements and intensities has indicated 
that no single method was always best. The C4 method usually gave best 
results when a system was intensifying or moving rapidly, but the FNWC 
prog, and often the low-order progs, appeared to give best results 
when the system was static or weakening. 

Investigation of the forecasts produced from the 17/OOOOZ and 
18/1200Z maps indicated that the above conclusions are applicable to 
other map times of this series. It does not follow, however, that 
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Figure 15. C4 24-hour 500-MB prog error. Contours are prog minus 
verifying analysis. Contour interval is 60 meters. 
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Figure 16. FNWC 24-hour 500-MB prog error. 
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these conclusions are applicable to all situations of the atmosphere.. 
High-index situations, for example, might be troublesome, although the 
high-order methods appeared to excel when forecasting large movement. 
Weak, summertime gradients might result in errors which are not 
apparent in this series There is little doubt, however, that the 
advection and conservation equations remain stable in typical meteoro- 
logical wind fields. Moreover, the theoretical phase and amplitude 
characteristics appear to remain valid for these fields 
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CHAPTER V 



SUMMARY 



Second-order and fourth-order advection equations have been 
developed for one dimension in conservation and advection forms. 
These equations have previously been shown to be stable for forward 



time steps when 



UAt 



. < 1.0, and in the case of the fourth-order 

Ax — 



conservation form, when 



f uAt 
' AX 



< 1.5. 



Phase retardation was also shown to be a minimum and amplitude 
damping a maximum for = 0.7. The advection form demonstrated 
somewhat better phase and amplitude characteristics than the 
conservation form. 

Use of Marchuk's method of fractional time steps permitted the 
one-dimensional equations to be applied to a two-dimensional finite- 
difference scheme without loss of stability or modification of phase 
and amplitude characteristics. 

In this study a controlled experiment utilizing a divergent, but 
linear velocity field was conducted. A known pattern representative 
of a simplified atmosphere was advected through 2ir radians by a 
velocity field of magnitude comparable to the atmosphere. Two hundred 
time steps were required. 

Positive, constant divergence was applied during the initial one 
hundred time steps; then negative divergence of equal magnitude was 
applied during the final one hundred steps. The object of the experi- 
ment was to determine how well the initial pattern was conserved. 
Crowley's theoretical stability analysis was substantiated by the 
experiment. 
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With a one-hour time step, computational instability of the 
Jacobian form destroyed the field before advection through it radians 
could be accomplished. Reduction of to one-half of the value 
used in the basic experiment and doubling the number of time steps 
eliminated the instability. It was found that the Jacobian form had 
better amplitude characteristics than the second-order advection 
method, but that the phase characteristics were inferior. 

The several advection methods were substituted for the advective 
term in a simple non-divergent barotropic model. The purpose of this 
experiment was to test the characteristics of the methods in a realistic, 
non-linear situation Operationally produced 500-MB analyses were 
obtained from FNWC for initial data and for forecast verification. 

Twenty- four and forty-eight hour forecasts of the 500-MB height 
field were produced using each form of the advective equation on eight 
consecutive map times beginning at 0000Z 15 January 1968. RMSE 
statistics were compiled for each of the methods on each of the map 
times. Machine plotted contour maps of all forecasts were generated 
for the map times 15/1200Z, 17/0000Z, and 18/1200Z Error charts were 
also plotted 

Analysis of the RMSE statistics indicated that the relative skill 
of the methods was Inversely proportional to the complexity of the 
method employed. Furthermore, all experimental methods except the 
fourth-order conservation form had lower RMSE than did the FNWC 
forecasts. 

Subjective analysis of the prognoses and error charts indicated 
that the inverse relationship was principally due to a tendency of the 
high-order equations to deepen existing low pressure centers excessively. 
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The poor performances of the FNWC progs was attributed to persistent 
errors in the western Pacific and eastern Asia. The cause of the FNWC 
errors is unknown. 

Detailed analysis of the 15/1200Z map time indicated that the 
fourth-order methods resulted in more rapid movement of troughs in 
strong gradient areas than occurred when lower order forms were used. 
Rapidly moving and deepening low centers were also handled better by 
the high-order methods, particularly the conservation method. Weakening 
low centers and static situations were generally handled better by the 
low-order methods. 

FNWC progs gave the best overall trough movement, primarily 
because of more rapid movement of the southerly portions. FNWC progs 
were at least equal to the lower order experimental methods when fore- 
casting for weakening cyclones. 

Crowley's equations were shown to be stable when subjected to the 
non-linear conditions of typical weather charts. The theoretical phase 
and amplitude characteristics appeared to remain valid under these 
conditions. 

No definitive conclusions concerning the relative skill of a 
baro tropic model employing the high-order advective equations could 
be made. A large number of trials under widely varying initial condi- 
tions would be necessary. Such further testing appears to be warranted 
and desirable in view of the favorable results of the limited experi- 
ments conducted thus far. 



58 



BIBLIOGRAPHY 



1. Crowley, Wo P„ "Numerical Advection Experiments Monthly 
Weather Review , XCVI (January, 1968), 1 - 11 

2. Haltlner, George J, "Numerical Weather Prediction,," Unpublished 
notes. Naval Postgraduate School, 1968, 

3. Leith, C, E, "Numerical Simulation of the Earth's Atmosphere," 
Methods in Computational Physics , Vol« IV, Academic Press, 1965 

4. Lewit, Howard L. "The Experimental Modification of a Barotropic 
Numerical Prediction Model," Unpublished Master's Thesis, Naval 
Postgraduate School, 1967„ 

5. Thompson, Phillip D, Numerical Weather Analysis and Prediction , 
MacMillan Company, 1961, 



59 



INITIAL DISTRIBUTION LIST 



No. Copies 



1. Defense Documentation Center 20 

Cameron Station 
Alexandria, Virginia 22314 

2. Library 2 

Naval Postgraduate School 
Monterey, California 93940 

3. Professor G. J. Haltiner 1 

Department of Meteorology 
Naval Postgraduate School 
Monterey, California 93940 

4. CDR. Dean R. Morford, USN 1 

2nd Weather Squadron 

Global Weather Central 

Offutt Air Force Base, Nebraska 68113 

5. Naval Weather Service Command 1 

Washington Navy Yard 
Washington, D. C. 20390 

6. Officer in Charge 1 



Navy Weather Research Facility 
Naval Air Station, Building R-48 
Norfolk, Virginia 23511 

7. Commanding Officer 1 

U. S. Fleet Weather Central 

COMNAVMARIANAS, Box 12 

FP0 San Francisco, California 96630 

8. Commanding Officer 1 

U. S. Fleet Weather Facility 
Box 72 

FP0 New York, New York 09510 

9. Commanding Officer 2 

Fleet Numerical Weather Central 
Naval Postgraduate School 
Monterey, California 93940 

10. Commanding Officer 1 

U. S. Fleet Weather Central 
Box 110 

FPO San Francisco, California 96610 



60 



1 



11. Commanding Officer 

U. S. Fleet Weather Central 
Box 31 

FPO New York, New York 09540 



12. Commanding Officer 1 

Fleet Weather Facility 
Navy Department 
Washington, D. C. 20390 

13. AFCRL - Research Library 1 

L. G. Hanscom Field 
Attn: Nancy Davis/Stop 29 
Bedford, Massachusetts 01730 

14. Superintendent 1 

Naval Academy 
Annapolis, Pferyland 21402 

15. Director, Naval Research Laboratory 1 

Attn: Tech. Services Info. Officer 
Washington, D. C. 20390 

16. Naval War College 1 

Newport, Rhode Island 02844 

17. Department of Meteorology 3 



Code 51 

Naval Postgraduate School 
Monterey, California 93940 

18. Department of Oceanography 1 

Code 58 

Naval Postgraduate School 
Monterey, California 93940 

19. American Meteorology Society 1 

45 Beacon Street 

Boston, Massachusetts 02128 

20. Bureau of Meteorology 1 

Attn: Library 
Box 1289 K, G. P. 0. 

Melbourne Vic. 3001 AUSTRALIA 

21. Program Director for Meteorology 1 

National Science Foundation 
Washington, D. C. 20550 

22. Office of Naval Research 1 

Department of the Navy 
Washington, D. C. 20360 



61 



23. Commander, Air Weather Service 2 

Military Airlift Command 
U. S. Air Force 

Scott Air Force Base, Illinois 62226 

24. Headquarters 2nd Weather Wing (MAC) 1 

United States Air Force 
APO New York, New York 09611 

25. Atmospheric Sciences Library 1 

Environmental Science Service Administration 
Silver Spring, Maryland 20910 

26. Director, Maury Center for Ocean Sciences 1 

Naval Research Laboratory 
Washington, D. C. 20390 



62 



JNCLASSiFIED 



Security Classification 



DOCUMENT CONTROL DATA - R & D 

( Security classification of title, body of abstract and indexing annotation must be entered when the overall report is classified) 



originating activity (Corporate author) 

Naval Postgraduate School 
Monterey, California 93940 



2a. REPORT SECURITY CLASSIFICATION 



UNCLASSIFIED 



2b. GROUP 



. REPORT T I TLE 

Numerical Experiments With High-Order Advective Equations 



J. DESCRIPTIVE NOTES (Type of report and m inclusive dates) 

Master's Thesis 



. AUTHOR(S) (First name, middle initial, last name) 

Dean R, Morford, Commander, USN 



. REPORT DA TE 

December 1968 



7a. TOTAL NO. OF PAGES 
60 



7b. NO. OF REFS 

5 



la. CONTRACT OR GRANT NO. 



b. PROJECT NO. 



9a. ORIGINATOR'S REPORT NUMBER(S) 



9fcfrO TH ER REPORT NO(S) (Any other numbers that may be as signed 
w this report 



0. DISTRIBUTION STATEMENT 



Distribution of this document is unlimited. 



11. SUPPLEMENTARY NOTES 



12. SPONSORING MILITARY ACTIVITY 

Naval Postgraduate School 
Monterey, California 93940 



3. ABSTRACT 



Second-order and fourth-order two-dimensional advective equations developed by 
W. P. Crowley were examined to determine their applicability to atmospheric models. 

One second-order and two fourth-order forms were evaluated from their perform- 
ances on a simple pattern advected by a linear, divergent velocity field. The same 
equations were substituted for the advection term of a simple barotropic forecast 
model to determine their performances on more general non-linear conditions. 

All forms of the equations remained stable in time and demonstrated the phase 
and amplitude characteristics predicted by Crowley The fourth-order "advection" 
form gave best results. 

When substituted in the barotropic model, the fourth-order forms lead to 
improved trough and low-center movements, but RMSE was slightly larger than that 
resulting from second-order forms The better RMSE of the lower order forms 
apparently resulted from their diffusive characteristics. 



DD 



FORM I47O 

1 NOV 65 I I KJ 

S/N 01 01 -807-681 1 



(PAGE 1 ) 



UNCLASSIFIED 



63 



Security Classification 



A-31408 



UNCLASSIFIED 



- Security Classification 



1 4 . 

KEY WO RDS 


L 1 N 


1 K A 


LINK B 


link c 


role 


W T 


ROLE 


WT 


ROLE 


W T 


Advection 

High-order equations 
Jacobian 

Numerical weather prediction 

. 












' 



DD , F .r.,1473 (BACK) UNCLASSIFIED 



S/N 0101-007-6821 



64 



Security Classification 



A- 3 t 409 



SH ELF BINDER 

Syrocuso, N. Y, 
: Sfocklon, Colif. 



